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Problem statement 


The dynamics and structure of crystal- melt interface are of great 
importance in crystal growth, melting, and other properties of 
solid liquid interface. The EAM , a semiempirical method developed 
by Daw and Baskes(l-3), for determining the energetics of atoms in 
a bulk enviroment , is a useful method for use in the computer 
simulation. The EAli functions of Pb in conjunction with the MD code 
are employed to study melting and disordering of Pb(100) surface. 
The EAM functions are fitted to the bulk experimental data at zero 
temperature and it is not clear if it will behave correctly at 
higher temperatures . In fact, it was shown(4) that use of a similar 
type of potential at higher temperatures might result in a negative 
coefficient of thermal expansion. The primary purpose of this 
report is to perform MD simulation combined with the EAM functions 
of lead to study melting of Pb(110) surface. This provide an 
accurate test of the EAM functions at higher temperatures . In 
particular, we calculate two dimensional structure factor, two 
dimensional pair distribution function, and density profile as a 
function of temperature up to the melting point of lead. 

Methodology 

Empirical pair potentials are suffering from two major problems, 
zero Cauchy pressure, i.e., P c /0 and equality of the vacancy 
formation energy E y with the cohesive energy E . In fact, in most 
of the metals P c =0 and E y <E c . To overcome these shortcomings the EAM 
functions were developed. In the EAM, the energy of an assembly of 
atoms is calculated as a function of their atomic coordinates. The 
major contribution to the energy of each atom comes from the 
embedding term F which can be interpretted as the energy that is 
gained when an atom is embedded into the charge density created by 
all the other atoms. Embedding function is supplemented by a two 
body potential $ mimicing the core-core repulsion, 

El=Fi+ - 5 | u *ij (r » ) ■ (l) 


where E: is the energy of atom i, F ^ is the embedding energy of atom 
i, p is the charge density at site i, is the electrostatic 
interaction between cores of atoms i and j, and r ^ is the 
interatomic distance between atoms i and j. p is approximated by 
the superposition of atomic charge densities of all the other sites 
except site i, 

Pr.S.Pj . ( 2 ) 

where is the atomic charge density. 
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F and $ ars determined by considering functional forms for 
them and fitting them to the bulk properties of solid. The 
exponential forms were considered for p and $(2). Parameters of the 
EAM functions of lead are tabulated in Ref. (3). Conditions under 
which simulations were performed are reported in table (1), where 
T is the temperature, and t^ and tj are equilibrium and total 
simulation times. 


In the MD simulation, one simply integrates the classical 
equations of motion from which velocity and position of the N 
particles in the system of volume V are calculated. Energy E 
remains constant in a microcanonical ensemble, so the time average 
of any property is equal to the time average over a microcanonical 
ensemble of configurations each having that same E, V, and 
|\j^0i'~cjQdieity)(5). Averages over a canonical ensemble, where V, T, 
and N are constant or over an isothermal - isobar ic ensemble, where 
T, P, and N are constant are more easily compared with the 
experimental results than averages over a microcanonical ensemble. 
Algorithms for T and P control have already been developed. 
Calculation of interparticle forces usually accounts for about 85% 
of the total simulation time. If all particles interact, 
calculation of the forces require N(N-l)/2 calls to the force 
routine. However, calls to the force routine can be reduced to N*N n 
if force is cut off beyond N n neighbor of each atom. 

We now briefly discuss the effects of hydrogen on metal -metal 
bonding using pair potential as well as the EAli and then compare 
our results. The interaction of the two metal atoms using the pair 
potential approach is, 

(3) 

where r is the interatomic distance. The force experienced by 
each atom is given by. 


We now put a hydrogen atom near to one of the metal atoms and 
calculate the potential and force between the two metal atoms; 


E 2 r H-t1' + ®H-H ^ r H-M ' ’ 

Therefore, in the pair potential approach, presence of the hydrogen 
does not change the force between the metal atoms. We now repeat 
the same calculation using the EAM formalism. 

Ep2F M ( p n ( r M _ M r M _ M .) , 

F 1 r 2F fi ( P N ) P N ( + r n-«) > 


(5) 

( 6 ) 


( 8 ) 


e 2 " f m ( P M (r M . H ) + p H (r H . M ) )+ F„( p M (r N . H ))+F H ( p H (r H . H )) + 

r H -H ) » ^ 


F 2" f h( P M +P H ) P«^ r h-M ^ + F h^ PH^M ( r H-^ + W r 
F 2 _F l =F t1^ Pf1 + PH> Ph ^ r l1-«^~ F «^ p fP Ph^h-^ ‘ 


H-M ^ ’ 


( 10 ) 
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Assuming Pu<<p M (which is a reasonable approximation), we can Taylor 
expand F 2 ~ F j> 

F 2“ F 1“ F M^ P|P Ph p M^ r IHp ’ ( - 12 - ) 

where F^>0, p H >0 , and pjj<0. This means that presence of the 
hydrogen will weaken the metal-metal bonding. 


Resu 1 ts 


a) Density profile p(z) 

The atom density per unit length is: 

p(2) = (l/^io>5 exp[-(z-z; ) 2 /2 o 2 ] 

where o- . 1 d no is the width parameter of the distribution, 
d 110 -a/2/2, and a-4.91 A. Result of the density profiles as a 
function of temperature is shown in Fig.(l). 

b) Structure factor Sj(T) 

The time-averaged squared modulus of the two dimensional 
structure factor for layer 1 at temperature T is: 

< | Sj| 2 >=<l/nj| S exp[-ig. r^ 2 >, 

where <> denotes the time averaging, n^ is the number of atoms in 
layer 1, and the sum is over the atoms of layer 1. <|Sj| > is 

calculated for g=2u/a(l ,0) . Our result for <|S 1 j / > as a function of 
temperature T for the top nine layers are shown in Fig. (2) 


Conclusions 

The maximum in p(z) for layer L2 is lower than for LI at 
temperatures 300 K and 400 K. A transfer of atoms between L2 to LI 
occurs between 400 K to 500 K. Above 500 K an adatom is formed 
above LI and leaves vacancies in LI and L2 . <[5j| > versus T curves 
show disordering of the top two layers begins near 400 K. The 
number of adatoms increases with the temperature and is maximum 
near the bulk melting temperature. 

Hydrogen embrittles the metal provided that the EAM approach 
is used for the energetics of the atoms. More work is needed to 
quantify the embrittlement near the grain boundaries in metals with 
and without the presence of hydrogen. 
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